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Abstract 

Background: About 6 million Americans suffer from heart failure and 70% of heart failure cases are caused by 
myocardial infarction (Ml). Following myocardial infarction, increased cytokines induce two major types of 
macrophages: classically activated macrophages which contribute to extracellular matrix destruction and 
alternatively activated macrophages which contribute to extracellular matrix construction. Though experimental 
results have shown the transitions between these two types of macrophages, little is known about the dynamic 
progression of macrophages activation. Therefore, the objective of this study is to analyze macrophage activation 
patterns post-MI. 

Results: We have collected experimental data from adult C57 mice and built a framework to represent the 
regulatory relationships among cytokines and macrophages. A set of differential equations were established to 
characterize the regulatory relationships for macrophage activation in the left ventricle post-MI based on the physical 
chemistry laws. We further validated the mathematical model by comparing our computational results with 
experimental results reported in the literature. By applying Lyaponuv stability analysis, the established mathematical 
model demonstrated global stability in homeostasis situation and bounded response to myocardial infarction. 

Conclusions: We have established and validated a mathematical model for macrophage activation post-MI. The 
stability analysis provided a possible strategy to intervene the balance of classically and alternatively activated 
macrophages in this study. The results will lay a strong foundation to understand the mechanisms of left 
ventricular remodelling post-MI. 



Background of myocytes results in significantly elevated interlukin-1 

Myocardial infarction is defined by pathology as myocytes (IL-1), tumour necrotic factor-/? (TNF-/J), IL-10, and 

necrosis and apoptosis due to prolonged ischemia. Since monocyte chemotactic protein- 1 (MCP-1) levels. MCP-1 

myocytes cannot divide and replace themselves, myocytes is a strong chemoattractants that recruit and confine 

in the infarct area deprived of oxygen die and are replaced monocytes to the injury site. It's been reported that over 

by a collagen scar. There is a series of cellular and molecu- 95% of monocytes differentiate to macrophages [2]. There 

lar activities respond to MI in the myocardium. Myocytes are two major types of macrophages post-MI: classically 

apoptosis appears in the first 6 to 8 hours post-MI, and activated macrophages (Ml) and alternatively activated 

necrosis occurs in 12 hrs to 4 days post-MI [1]. Necrosis macrophages (M2). Porcheray has reported a switch 

between Ml and M2 macrophages with in vitro stimuli 
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transition [3]. Since Ml and M2 macrophages are 
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responsible for extracellular matrix (ECM) destruction and 
construction, respectively, the transition and dynamic bal- 
ance between two macrophage phenotypes might lead to 
the balance between ECM destruction and construction, 
and thus determine the ECM remodeling post-MI.[4] 
Therefore, characterizing macrophage activation pattern 
is essential to better understand the ECM remodeling 
post-MI. 

A large amount of experimental research has been con- 
ducted to elucidate the underlying mechanisms of 
macrophage activation, and an abundant accumulation of 
experimental results define on macrophage responses to 
different stimuli. There is a need, however, to systemi- 
cally analyze the accumulated data and integrate the 
results into a framework that will allow a more complete 
understanding. To address this need, several mathemati- 
cal models have been established to characterize the 
effects of macrophages on wound healing, inflammatory 
responses, and collagen synthesis post-MI [5-9]. How- 
ever, most models do not consider the effect of macro- 
phage activation patterns and ignore the differences 
between macrophage phenotypes. Therefore, the aim of 
this study was to establish and validate a set of ordinary 
differential equations to characterize macrophage activa- 
tion patterns post-MI. Since our mathematical model 
was established based on in vivo and in vitro experimen- 
tal results, all parameters in the model were determined 
by the averages of the experimental data. 

Results 

We have collected experimental data from adult C57 mice 
and built a framework to represent the regulatory relation- 
ship among cytokines and macrophages. Based on this fra- 
mework, we established a set of nonlinear differential 
equations to characterize the regulatory relationship for 
macrophage activation in the left ventricle post-myocardial 
infarction using physical chemistry laws. Our framework 
and the mathematical model were established based on 
the following three assumptions. 

1) All monocytes that migrate to the infarct region 
are differentiated to unactivated macrophages [10]. 

2) All activated macrophages are differentiated from 
unactivated macrophage since previous studies have 
shown that <5% of macrophages undergo mitotic 
division [11,12]. 

3) All parameters and coefficients in this model are 
constant. 

Framework of regulatory relationship for macrophage 
activation 

In this framework, myocytes and monocytes were con- 
sidered as inputs to the system. Cellular densities of Ml 
and M2 macrophages were considered as the outputs of 



the system. We chose IL-1, IL-10, and TNF-a as three 
molecules which regulate macrophage activations in this 
mathematical model since they were well-recognized as 
stimuli for macrophage activation [13-15]. Ml macro- 
phages and myocytes secrete IL-1 and TNF-a. M2 
macrophages secrete IL-10 [3]. Further, TNF-a and IL-1 
promote Ml activation and IL_10 promotes M2 activa- 
tion [16]. IL-10 inhibits TNF-a, IL-1, and itself [17]. 
The input-output and regulating relationship were 
shown in Figure 1. 

Input to the framework 

Temporal profiles of monocytes and myocytes densities 
were used as inputs to our mathematical framework 
(Figure 1). Myocytes density in healthy adult mice was 
6 x 10 9 cells/ml as an initial value. Myocytes numbers 
monotonically decreases post-MI and is directly asso- 
ciated with LV wall thickness. We have measured the 
LV wall thickness at days 0, 1, 3, 5, and 7 post-MI. The 
temporal profile of myocytes was determined by com- 
bining the initial value and the monotone progression 
trend (the crosses) as shown in Figure 2(A). 

Macrophages density in the left ventricle of healthy 
adult mice is 2000 cells/ml, which will be used as initial 
values of unactivated macrophage density in this study 
[18,19]. Yang et al has measured temporal profiles of 
macrophages post-MI in mice at days 1, 2, 4, 7, 14, 21, 
and 28 [20]. The temporal profile of macrophages was 
obtained by fitting the experimental results as a continu- 
ous function as shown in Figure 2(B). 

Based on our assumptions 1 and 2, all macrophages 
were differentiated from monocytes and emigrated from 
infarct area to the lymph node system. The estimation 
of unactivated macrophage based on the experimental 
results [18] is shown as follows, 

M un {t) = M un {0)e-^ + f e-^Mds, (1) 
Jo 

where M denotes the differentiation rates of mono- 
cytes, M un denotes the unactivated macrophage density, 
and ft denotes the emigration rate of inactivated macro- 
phage. Based on the temporal profile of unactivated 
macrophage, the monocytes differentiation rate can be 
estimated as shown in Figure 2(C). 

Mathematical model for macrophage activation 

The mathematical model of macrophage activations is a 
set of nonlinear differential equations represented by 
cellular densities (cell number/ml) of M u „, Ml and M2, 
and concentrations (pg/ml) of chemical factors such as 
IL-1, IL-10, and TNF-a. Cellular densities were deter- 
mined by the difference between immigration and emi- 
gration rates. Concentrations of chemical factors were 
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Figure 1 A framework representing regulatory relationship macrophage activation post-MI. Myocytes and monocytes are two inputs as 
shown in the blue dotted box. The red dotted box denotes outputs of the system, including M1 and M2 densities. Secretion function is 
denoted by thick arrow, activation and differentiation are denoted by thin arrow. The regulatory relationship shown in the framework was 
determined by published results [3,13-17]. 



determined by the balance between their synthesis and 
degradation rates. The established mathematical model 
can be written as 
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where M un , My M 2 denote the cell densities of unacti- 
vated macrophages, Ml macrophages, and M2 macro- 
phages, respectively. Variables IL W , T w and IL 1 denote 
the concentrations of IL-10, TNF-a, and IL-1. Variable 
denotes M the differentiation rate of monocytes and M c 
denotes the myocytes density. The parameters used in 
these equations with their biological meanings, experi- 
mental values, units, and references are listed in Addi- 
tional file 1. All parameters were determined based on 
the published data or estimation from other mathemati- 
cal models [3,6,21-28]. 
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Figure 2 Temporal profiles of monocytes and myocytes. The wall thickness of left ventricle and observed a monotone decrease of wall 
thickness at day 1, 3, 5, 7 post-MI, suggesting the same trend of myocytes density based on the experimental results [18]. Temporal profiles of 
myocytes and the monotone progression trend are shown in Figure 2(A). We established the temporal profiles of monocytes based on our own 
experimental measurement on total macrophage density and the control value of macrophage density in health adult mice is 5 x 1 0 3 cells/ml. 
Temporal profiles of monocytes are shown in Figure 2(B). The estimation of monocytes for each day are shown in Figure 2(C). 



Equation 2 determines the density of unactivated 
macrophages in the infarct area. For the construction 
part, the unactivated macrophages are differentiated 
from monocytes as shown in Figure 2(C). For the 
destruction part, the unactivated macrophages are acti- 
vated to Mi or M 2 - Additionally, inactivated macro- 
phages do not die locally in the scar tissue but die out 
in the lymph node system [27]. 

Equation 3 determines the activation rate of Ml macro- 
phages. For the construction part, IL-1 and TNF-a pro- 
mote Ml activation. Parameters k 2 and k 3 denote the 
activation rates of Ml macrophages by IL-1 and TNF-a. 
Hill equations are used to represent the promotion effects 
of IL-1 and TNF-a and parameters c IL \ and Cj„ are the 
effectiveness of IL-1 and TNF-a promotion on Ml calcu- 
lated based on the experimental results [3,29] . Steinmuller 
has shown the transition between Ml and M2 phenotypes 



in vivo [21]. Correspondingly, we use parameter k\ to 
denote the transition from Ml to M2 and parameter k\ for 
the transition from M2 to Ml [21]. The destruction part 
includes emigration of macrophage (^)and transition from 
Ml to M2 macrophages (ki) [27]. 

Equation 4 determines the activation rate of M2 
macrophages. The construction part is denoted by acti- 
vation of M2 macrophages promoted by IL-10, and 
transition from Ml to M2. IL-10 promotes M2 activa- 
tion and this activation rate has been approximated by 
parameters k 4 based on the experimental results [3]. 
The transition rate from Ml to M2 is denoted as k\. 
The destruction part includes emigration of M2 macro- 
phages (fi) and transition from M2 to Ml macrophages 
(fe'j), similarly as described in equation 3. 

Equation 5 determines the secretion rate of IL-10. For 
construction part, IL-10 is secreted by M2 macrophages, 
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and parameter k s denotes the secretion rate of IL-10 by 
M2 macrophages [22,30]. The destruction part includes 
the self-inhibition and degradation of IL-10. A Hill 
equation is employed to represent this self-inhibition 
effect and parameter c± denotes the self-inhibition effect 
of IL-10 post-MI [6,17]. Parameter dn m denotes the 
decay rate of IL-10 determined by its half-life time [13]. 

Equation 6 determines the deposition rate of TNF-a, 
which is secreted by both Ml macrophages and myo- 
cytes [6,23,28,31]. We used in vitro results from Meng 
to determine the secretion rate of TNF-a by Ml and 
secretion rate of TNF-a by myocytes (A) is determined 
by Horio's experimental results [23,31]. The inhibition 
of IL-10 is presented by a Hill equation where para- 
meter c represents the effectiveness of IL-10 inhibiting 
TNF-a [25]. The destruction part is denoted by the 
degradation of TNF-a. Parameter dj a is the decay rate 
of TNF-a determined by its half-life time [15,26]. 

Equation 7 determines the deposition rate of IL-1. IL-1 
is secreted by both Ml macrophage and myocytes. Para- 
meter k 7 denotes the secretion rate of IL-1 in cultured rat 
cardiac myocytes [31]. The inhibition of IL-10 is pre- 
sented by a Hill equation similarly as in equation 6 [25] . 
In the destruction part, parameter dui represents the 
decay rate of IL-1 determined by its half-life time [14]. 

Stability analysis 

If there is no myocardial infarction, monocytes differen- 
tiation and myocytes apoptosis should be at a very low 
level, and the studied macrophage activation pathway 
should maintain homeostasis. We have calculated the 
equilibrium point of the system without any input and 
performed Lyapunov stability analysis. Our analysis 
showed that without any monocytes differentiation and 
myocytes secretion, the system would stay at the origin 
when t — > oo. 

In the case of myocardial infarction, myocytes apopto- 
sis and necrosis triggered inflammatory responses and 
significant monocytes differentiation, which will drive the 
system to a new equilibrium point. Correspondingly, the 
cell densities of Ml and M2 increase post-MI. We have 
obtained a steady state as E = [20, 1200, 3500, 0.73, 1.1, 
5.9] from our computational simulations. The steady 
states match with the experimental measurements col- 
lected from healthy adults without myocardial infarction 
[32]. In addition, the stable region of the established 
mathematical model depends on the strength of the 
input. 

Computational results 

Computational simulations of macrophage activation 
were carried out by solving the nonlinear differential 
equations with MATLAB. The initial conditions of 
unactivated, Ml and M2 macrophage densities were 



chosen as M un = 2000 cells/ml and Mi(0) = M 2 (0) = 0 
cells/ml. The concentrations of IL-1, TNF-a, and IL-10 
were set as 0.1 pg/ml. The inputs of this system were 
shown in Figure 2. Outputs of the system, Ml and M2 
densities, were shown in Figure 3. The concentrations of 
IL-1, IL-10 and TNF-a were shown in Figure 4. Our 
computational results were shown as solid lines while 
the experimental results were shown as discrete crosses 
in these figures [17,33,34]. 

Our computational results demonstrated that from 
days 0 to day 3 post-MI, cellular densities of the Ml 
phenotype increased at a faster rate than the M2 pheno- 
type. At day 10, the M2 phenotype dominates over the 
Ml phenotype. This prediction agrees with the results 
reported by Troidl [4]. Additionally, temporal profiles of 
IL-1 and TNF-a significantly increased from days 0 to 1 
post-MI in our computational simulations, which match 
the experimental results reported by Sumitra [33]. Com- 
parison between the computational and experimental 
results demonstrates a similar trend for the temporal 
profiles, suggesting the effectiveness of our model. 

Discussion 

This study established a mathematical model for macro- 
phage activation in the left ventricle post-MI by combin- 
ing experimental and computational approaches. This is 
the first mathematical model focusing on the dynamic 
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Figure 3 Temporal profiles of cell densities of Ml (red line), 
M2 (blue line), and the total macrophage (black line) in LV 28 
days post-MI. Computational results were normalized to initial 
conditions and shown in solid lines. The initial cells number of 
unactivated macrophage is 2 x 1 0 3 cells/ml. The differentiated rate 
from monocytes to M1 macrophages a is 0.5. Then, the initial 
condition of macrophage densities were chosen as /M1(0) = M2(0) = 
0 cells/ml. Previously published experimental results were 
normalized to the corresponding measurements in control group 
and are shown as black crosses collected from C57 adult mice at 
days 1, 3, 5, 7, 14, 21, 28 [35]. Values are mean ± SE. All experiments 
were carried out in mice with Ml. 
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Figure 4 Concentrations of IL-10, IL-1 and TNF-a from day 0 to day 30 post-MI. The comparisons between the experimental 
measurements and temporal files were presented. In Figure 4 (A), experimental measurements of IL-10 from adult C57 mice at day 0.25, 3, 5 
post-MI [34], were shown as crosses with standard deviations shown as the bar. In Figure 4 (B), experimental measurements of IL-1 from adult 
mice at day 0.25, 1, 3, 5, 7, 21 [36], were shown as crosses with standard deviations shown as the bar. In Figure 4 (C), experimental 
measurements of TNF-a from adult rats at day 0.25, 1, 3, 5, 7, 14, 21 [37], were shown as crosses with standard deviations shown as the bar. 



interactions among cytokines, myocytes, monocytes, and 
macrophages. Computational predictions based on this 
mathematical model match with experimental measure- 
ments, suggesting effectiveness of the model. In addition, 
our stability analysis provided insight for the activation 
pattern of macrophages post-myocardial infarction. In 
our mathematical model, there are two inputs, myocytes 
apoptosis and monocytes differentiation. In this study, 
we have predicted a stable equilibrium for homeostasis, 
which means without myocardial infarction or following 
small injury stimuli, macrophage densities and concen- 
trations of IL-1, IL-10, and TNF-a should stay at the 
equilibrium. After MI, the monocytes differentiated into 
macrophage and apoptotic myocytes secreted significant 
amounts of cytokines to activate the macrophages. The 
strength of the monocytes differentiation and myocytes 
apoptosis (inputs of the system) drive the system to 



different states while all states will be bounded due to the 
bounded strength of the inputs. 

However, there exist some differences between compu- 
tational predictions and real experimental results. To 
address this issue and the variation in different experi- 
ments, a stochastic parameter distribution will need to be 
introduced to replace the constant parameters. In addi- 
tion, more detailed measurements on monocytes and 
myocytes and concentration temporal profiles of IL-1, 
IL-10 and TNF-a from C57 mice will be carried out in 
our future research, which will help to solidify our 
equations. 

Conclusions 

Our study has established framework for macrophage 
activation and used ordinary differential equations to 
model the cellular interactions between macrophage 
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activation types post-MI. The results on stability analysis 
can be used as a useful tool to predict the behaviour of 
biological systems. 

Methods 

To incorporate the experimental data, curve fitting algo- 
rithm was applied to obtain temporal continuous density 
profiles of myocytes and monocytes based on discrete 
experimental data. 

Stability analysis of the established mathematical model 

To analyse the stability of the proposed mathematical 
model, we have calculated the equilibrium point of the 
system and performed the Lyapunov stability analysis of 
the system. 

In our mathematical model, equations (2-7) are 
six first-order equations with input M and M c . 
Now we use x x , x 2 , x 3 , x 4 , x 5 , x 6 to denote M un , 
Mi, M 2 , IL 10 , T a , IL\. We first examined the stability 
of the system without any input and obtained an equi- 
librium of Ei = {x\, x* 2 , x\, x\, x* 5 , x* 6 ) = (0,0,0,0,0, 0). 
With the temporal input of monocytes differe- 
ntiation and myocytes secretion, our computa- 
tional simulation generated a steady state as 
E 2 = (x*,x* 2 ,xl,x* 4 ,x* 5 ,xl) = (20, 1200,3500,0.73, 1.1,5.9). To 
further analyze the stability property of the system, we 
chose a positive definite Lyapunov function 

and obtained its derivative as the following equation 

V(x) = 2_, j x i x i 

= —ax] - X\ \k 2 x< * 6 + k,xi— — + k 4 x 1 —^ i — 1 

- (n + ki)x\ + k 1 x 2 xij^- i + k 3 xix 2J ^- + k\x 3 x 2 (9) 

- (fe'i + /j.)xj + hx 2 x 3 + heix 3I ^^ 

+ k5X3X4 7?ixl ~ d ^ X 5 + k <> X 2 X 57^ - dl Ll xl + k 7 X 2 X 6 ^. 

By applying the boundary of the Hill equations, we got 

V(x) < —(// + ^2+^3+ fe 4 )x^ — (// + k\)x\ + k2*2*l + ^3^1X2 + ^1X3X2 

- (fe'i + ti)x\ + k 1 x 2 x i + fe 4 xix 3 - d, Lm xl + fe 5 X3X 4 - ii To x| + fe 6 X 2 X5 (10) 

- d, L1 xl +k 7 x 2 x 6 . 

Applying the parameters in Table 1 to equation (10), 
we got 

V(x) < -1.6xJ - 0.275l| - 0.25*! ~ 2 5x 4 ~ 55x 5 ~ 10 5x 6 + l lxiU + 0.3x,x 3 
0.125x 2 x 3 + 0.0005x 3 x 4 + 0.0007X2*5 + 0.0006x 2 x 6 

< -0.02xJ - 0.02;t! - 0.05x1 - 0.5x| - 5x\ - 0.5x\ - 1.28(xi - 0.429x 2 ) 2 /1 -1 \ 

- 0.3(x, - 0.5x 3 ) 2 - 0.25(x 3 - 0.25x 2 ) 2 - 2(x„ - 0.000125x 3 ) 2 

- 50(x 5 - 0.000007x 2 ) 2 - 10(x 6 - 0.00003x 2 ) 2 

< 0. 

Since the derivative of Lyapunov function is negative 
semi-definite and the semi-definite is satisfied with all 
states equal to zero, the system is globally asymptotically 



stable without any input. Given any bounded differenti- 
ate rate of monocytes and myocytes density as input, 
the system will have bounded states. 

Additional material 



Additional file 1: Table 1. Pre-determined parameters from literature 
search. 
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